## Measuring Political Support and Issue Ownership Using Endorsement Experiments,
## with Application to Militant Groups in Pakistan

library(R2jags)

odd <- function(x) x%%2==1

J <- 4 ## number of issue areas asked about
K <- 5 ## number of endorsement groups (including control)
n.province <- 4
n.division <- 26
province <- c(1,1,1,1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,3,3,4,4,4,4,4,4)
n.per.division <- c(118,0,313,403,579,495,208,131,203,473,311,0,293,0,84,287,50,0,215,288,103,0,210,320,67,61) 
N <- sum(n.per.division) ## number of survey respondents
n.age <- 4

  ## Ideal Points
  delta.division <- c(-1,-1,-0.75,-0.75,-0.5,-0.5,-0.25,-0.25, #province 1
                  -0.5,-0.25,0,0.25,0.5,  #province 2
                  -0.75,-0.5,-0.25,0,0.25,0.5,0.75,  #province 3
                  0,0,0.5,0.5,1,1)  #province 4
  division <- NULL
  for (i in 1:n.division){  division <- c(division, rep(i,n.per.division[i])) }
  
  delta.age <- c(0.5, 0.5, 0, -1)
  age <- sample(1:4, size=N, replace=T)
  
  delta.female <- 0.5
  female <- rbinom(N, 1, p=0.47)
    
  x <- rep(NA, N)
  for (i in 1:N){ x[i] <- rnorm(1, mean=(delta.age[age[i]] + delta.female*female[i] + delta.division[division[i]]), sd=1) }
  x <- x / sd(x)
  
  beta <- runif(J, 0.5, 2) ## discrimination parameter
  alpha <- matrix(0.1,J,4)
  for (j in 1:J){  ## variation in thresholds, but distance ensured
	  while(alpha[j,2]-alpha[j,1] < 0.5 | alpha[j,3]-alpha[j,2] < 0.5 | alpha[j,4]-alpha[j,3] < 0.5){
		  tmp <- rnorm(4, mean=0, sd=1)
		  alpha[j,] <- tmp[order(tmp)]
	  }}

  lambda.division <- matrix(c(rep(0,n.division), rep(0, n.division), delta.division, -delta.division, delta.division/3), n.division, K)
  theta.age <- matrix(c(rep(0,n.age), 0.5,0.25,0,-0.75, 0.25,0,0,-0.25, -0.75,-0.25,0.25,0.75, 0,0,0,0), n.age, K) ## for each K, values must sum to 0
  theta.female <- c(0, 0.25, 0.5, 0, -0.25)
  
  omega <- c(0, 0.3, 0.3, 0.5, 0.5)
  
  s <- array(NA, c(N,J,K)) ##Individual level impact of group endorsement k
  for (i in 1:N){
    for (j in 1:J){
      for (k in 1:K){
        s[i,j,k] <- rnorm(1, mean=(theta.age[age[i],k] + theta.female[k]*female[i] + lambda.division[division[i],k]), sd=omega[k])
  }}}

  
##Ordered Responses (5 choices)
Y <- array(NA, c(N,J,K))
for (i in 1:N){
	treat <- ifelse(odd(i),TRUE,FALSE) ## Half of Jake's respondents receive the control, remaining half split among treatments
	for (j in 1:J){
		for (k in 1:K){
			p1 <- pnorm(alpha[j,1] - beta[j]*(x[i] + s[i,j,k]))
			p2 <- pnorm(alpha[j,2] - beta[j]*(x[i] + s[i,j,k])) - (p1)
			p3 <- pnorm(alpha[j,3] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2)
			p4 <- pnorm(alpha[j,4] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2 + p3)
			p5 <- 1 - (p1 + p2 + p3 + p4)
			
			p <- c(p1,p2,p3,p4,p5)
			Y[i,j,k] <- sample(c(1:5), size=1, prob=p) ## Full array of potential outcomes
		}}
	if (treat==FALSE) { Y[i,,2:K] <- NA}
	if (treat==TRUE) {
		Y[i,,1] <- NA 
		for (j in 1:J){ 
			answered <- sample(2:K,1); Y[i,j,-answered] <- NA  
		}
	}}

Y2 <- matrix(NA,N,J)
endorser <- matrix(NA, N,J)
for (i in 1:N){
	for (j in 1:J){
		endorser[i,j] <- (1:5)[!is.na(Y[i,j,])]
		Y2[i,j] <- Y[i,j,endorser[i,j]]
	}}
Y <- Y2


  lambda.division.real <- lambda.division
  theta.age.real <- theta.age
  theta.female.real <- theta.female
 
  #######################################
  ##  Prep and send data through JAGS  ##
  #######################################
  data <- list ("N", "J", "K", "Y", "endorser",
    "division", "n.division", 
    "n.age", "age", "female",
    "province", "n.province",
    "lambda.division.real", "theta.age.real", "theta.female.real") 
  inits <- function() {list(alpha0=matrix(seq(-2,2,length.out=(4*J)),4,J), beta=runif(J, 0.2, 2), 
    x=rnorm(N), s=matrix(rep(0, N*J), N,J), 
    delta.division=rep(0,n.division), delta.province=rep(0,n.province),
    delta.age.raw=rnorm(n.age), delta.female=rnorm(1),
    xi.age=runif(1),
    theta.age.raw=matrix(c(rep(0,n.age), rnorm(n.age*(K-1), n.age, K)), n.age, K), theta.female=c(0,rnorm(K-1)),
    lambda.division=matrix(c(rep(0,n.division), rep(0,n.division*(K-1))), n.division, K), theta.province=matrix(c(rep(0,n.province), rep(0,n.province*(K-1))), n.province, K),
    xi2.age=c(0,runif(4)),
    omega=c(0,runif(K-1)) )} 
  params <- c("lambda.division.error", "theta.age.error", "theta.female.error") 
  model <- "ModelSimulationsDivisionWithICovariates.txt" 
  fit <- jags(data, inits, params, n.iter=20000, model.file=model)
  fit <- autojags(fit, n.iter=20000, n.thin=20, Rhat=1.02, n.update=5)
  save.image(file = Outfile)
  
  